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ABSTRACT 

We present a new fast and efficient approach to model structure formation with aug- 
mented Lagrangian perturbation theory (ALPT). Our method is based on splitting the dis- 
placement field into a long and a short range component. The long range component is com- 
puted by second order LPT (2LPT). This approximation contains a tidal nonlocal and nonlin- 
ear term. Unfortunately, 2LPT fails on small scales due to severe shell crossing and a crude 
quadratic behaviour in the low density regime. The spherical collapse (SC) approximation has 
been recently reported to correct for both effects by adding an ideal collapse truncation. How- 
ever, this approach fails to reproduce the structures on large scales where it is significantly less 
correlated with the A^-body result than 2LPT or linear LPT (the Zeldovich approximation). 
We propose to combine both approximations using for the short range displacement field the 
SC solution. A Gaussian filter with a smoothing radius rs is used to separate between both 
regimes. We use the result of 25 dark matter only A^-body simulations to benchmark at z = 
the different approximations: 1st, 2nd, 3rd order LPT, SC and our novel combined 2LPT-SC 
model. This comparison demonstrates that our method improves previous approximations at 
all scales showing ~25% and ^75% higher correlation than 2LPT with the A'-body solution 
at /c = 1 and 2 h Mpc~^, respectively. We conduct a parameter study to determine the optimal 
range of smoothing radii and find that the maximum correlation is achieved with rs = 4 — 5 

Mpc. This structure formation approach could be used for various purposes, such as 
setting-up initial conditions for A^-body simulations, generating mock galaxy catalogues, per- 
form cosmic web analysis or for reconstructions of the primordial density fluctuations. 
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1 INTRODUCTION 

Numerical simulations are now an indispensable tool to study cos- 
mological structure formation. The Poisson-Vlasov equations, de- 
scribing the gravitational interaction of matter in an expanding 
Universe, are usually solved with a large number of test parti- 
cles using_so^ned_iV2^^ reviews on this subject 
see iBertschingedll 9981 : iKuhlen et al 201 2|). The m ost extended nu- 
merical techniques are the tree codes " fearnes & Hut 1986.) and 
the adaptive particle-mesh (PM) methods jKravtsov et al. 1 19971 : 
iTev ssier '2002). An example for a hybrid tree-PM code is GAD- 
GET ( Springel 2005) which uses the PM method to evaluate long 
range forces and the tree method for short range interactions. The 
forces between particles need to be recomputed in each time inter- 
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val, requiring in total thousands of time steps to complete a single 
simulation. 

Lagrangian perturbation theory (LPT) gives an approximate 
solution to structure formation by computing the displacement 
field from the primordial density fluctuations to the present Uni- 
verse in one step (for a review see Bernardeau et al. 2002). It 
has a wide range of applications in cosmology, giving an analyt- 
ical understanding to structure formation (.Zel ' dovich. J 97 0.) . Cur- 
rently, it is commonly used to set-up initial condi t ions for numer - 
ical A^-body simulations (see ICrocce et al.ll2006l : I Jenkins '2010*). 
Another field of appl ications is the reconstruction of peculiar ve- 
locity fields (see e.g . iMohavaee & Tullvll2005l : iLavaux et al.ll2008l : 
iK itaura et al. 1201 2|) and primordial density fl u ctuations (see e.g. 



Moha yaee et all l2006l : lEisenstein etaP l200l iKitaura & Angulol 
120121 : lKitau^ l2012'). It can also be used to change the cosmol- 
ogy of A^-body simulations jAn gulo & Whit9 r2010). The need for 
simulations to estimate uncertainties in cosmological parameter es- 
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Figure 1. Cell-to-cell correlation between the linear initial overdensity field Di6(q, z = 100) and the corresponding approximations for the divergence of 
the displacement field for the 10th realisation of our set of simulations. The solid black line represents the LPT/Zeldovich approximation and the green curve 
the local SC model, which approximately fits the mean A^-body relation. The nonlocal relations are given by the contours for various approximations: left 
panel: 2LPT (quadratic relation), middle panel: 3LPT (cubic relation) and right panel: combined 2LPT-SC with rs = 4 Mpc//i. The dark colour-code 
indicates a high number and the light colour-code a low number of cells. 



timations from large scale structure surveys has raised the inter- 
est in approximate effici ent structure formation models, which can 
be ma ss ively used. See Scoccimarro & Shethl (l2002h : lMonaco et al.l 
(l2002h : iManera et al. (2012 ) for generation of mock galaxy cata- 
logues: and Schneider et "aLl(l201 Ih to increase the volume of a set 



of A^-body simulations with smaller volumes. An improvement to 
linear LPT is given by second order LPT (2LPT) including non- 
local tidal field corrections. However, this is known to be a poor es- 
timator in hig h and low density regions, being strongly limited by 
shell crossing (ISahni & Shandarinll996l : lNevrinckl2013h . Recently, 
local fits based on the spherical collapse model have been proposed 
dMohayaee et al. 2006; Neyrinck 2013), which better match the 
mean stretching parameter (divergence of the displacement field) 
of A/^-body simulations. We propose in this work to combine the 
superiority of 2LPT on large scales with the more accurate treat- 
ment of the spherical collapse (SC) model on small scales includ- 
ing a collapse truncation of the stretching parameter, which acts as 
a viscosity term. Our algorithm splits the displacement field into a 
long-range and a short-range component, the first one being given 
by 2LPT and the second one by the truncated SC model. Both are 
combined by using a Gaussian filter with smoothing scales of 4-5 
Mpc/h radii, being this scale the only free parameter in our model. 



2 THEORY 

Let us define the positions of a set of test particles at an initial time 
ti by q and call them the Lagrangian positions. The final comoving 
positions x (called Eulerian positions) corresponding to the same 
set of test particles at a final time t / are related to the Lagrangian 
positions q through the displacement field, ^: 



(1) 



Hence, the displacement field encodes the whole action of gravity 
during cosmic evolution. An approximation is to consider that the 
displacement field is a function of the initial conditions only, and 
can be described by straight paths. The various models consider 
the relation between the divergence of the displacement field and 
the linear initial field: = '0(^*^^^) = V • *(^*^^^), where ^ is 
the so-called stretching parameter. Let us call the previous equation 
the stretching parameter relation. Lagrangian perturbation theory to 



third order yields the following expression for cu rl-free fields (see 
lBuchertlll994l : iBouchet et al.ll 19951 : ICatelanlll995h : 



V^SLPT = V • ^3LPT 



(2) 



where Di is the linear growth factor, D2 the second order growth 
factor, {Dsa, i^3b, ^3c} are the 3rd order growth factors corre- 
sponding to the gradient of two s calar potentials (0a , 0k ). Partic- 



ular expressions can be found in lBouchet et alj (Il995h and lCatelanI 
(1995): D2 = Dsa = -1/3 

Dsb = 1/4 • 10/21 6^^\q) represents 

the 'second-order overdensity' and is related to the linear overden- 
sity field by the following quadratic expression: 



<5<^)(q) = ^(.^«(q)<^«(q)-[.^«(q)]^), 



(3) 



i>j 



The potentials (f)^ ^ and (f)^ ^ are obtained by solving a pair of Pois- 
son equations: Vg(/)*^"^'*(q) — d^-^^q), where d^-^^q) is the linear 
overdensity, and Vg^*^^-* (q) = 6^'^\q). The first term is cubic in 
the linear potential 



(4) 



and the second term is the interaction term between the first- and 
the second-order potentials: 



^(3) 



(5) 



(see lBucherllll994l : iBouchet et al1ll995l: ICatelanl[l99l . Keeping 
terms only to firs t order is called the Zeldovich approximation 
(IZerdovich|[T97Qh and keeping terms to second order yields the 
2LPT approximation. 

Based on the nonlinear s p herical collapse approx imation, de- 
veloped bv lBemardeaul (ll994l) . lMohavaee et"aD(l2006h found a lo- 
cal nonlinear expression for the stretching parameter relation: 



1/2 



1 



(6) 



This analytic formula has been recently found to fit very well the 
mean stretching parameter relation from an A/^-body simulation 
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Figure 2. Slices for various structure formation models atz = Oof~13/i~^ Mpc thickness (averaged over 9 neighbouring slices) through a volume 
of 180 Mpc per side showing the logarithm of the density field (ln(2 + 5)) using the same colour range ([—1, 1.5]) on a mesh with 128^ cells after 
applying clouds-in-cells. Upper left panel: 2LPT; upper right panel: A^-body; lower left panel: combined 2LPT with 2LPT including collapse threshold 
(2LPT-2LPTC); lower right panel: combined 2LPT with SC model (2LPT-SC) for the 10th realisation of our set of simulations. 



(lNevrinckll2013l) . We will call t he spherica l colla pse (SC) approx- 
imation the improved model bv lNevrinckl (EoBi) . which includes 
an ideal collapse limit. The stretching parameter cannot be smaller 
than —3, considering that for the ideal collapse to a point we have 
V • X = and thus = V • * = -3. 



3 METHOD 

The method we present here is straightforward and computation- 
ally extremely fast. Our ansatz is to expand the displacement field 
^(g, z) into a long range ^L(g, z) and a short range component 

^[ci,z) = *L(q,2;) + *s(g,2:). (7) 

We will rely on 2LPT for the long range component: ^A(g, z) — 
^2LPT(g, z). It is known that 2LPT fails towards small scales (see 
next section). Therefore, we propose to filter the displacement field 



resulting from the 2LPT approximation. We consider in this work 
a Gaussian filter /C(g, rs)= exp ( — |g|^/(2rs)), with rs being the 
smoothing radius. 

*L(qr,2;) = /C(g,rs)o*A(g,2;). (8) 

The spherical collapse approximation yields an extremely good lo- 
cal fit to the mean relation between initial and final stretching pa- 
rameter. However, as it does not take into account the nonlocal tidal 
component it dramatically fails on large scales. We will use this 
approximation to model the short range component (^3(9, ^) = 
*sc(<7,^)): 

*s(g, z) = *B(g, z) - /C(g, rs) o *b(<7, z) . (9) 

One may improve the method finding better approximations for the 
long range or the short range components. We dub the 
computer-code based on our method: augmented Lagrangian per- 
turbation theory ALPT. 



© 0000 RAS, MNRAS 000, 000-000 



4 F. S. Kitaura and Steffen Heft 



4 NUMERICAL EXPERIMENTS 

We benchmark different structure formation approximations by 
comparing them with the results from A/^-body simulations taking 
25 randomly seeded initial conditions. The simulations were done 
with GADGET- 3, an improved version of the publicly available cos- 
mological code GADGET-2 (last described in Springel 2005), using 
384^ particles on a box with 180 Mpc side. The initia l condi- 
tions were generated using the WMAP7 parameters (K omatsu et al.l 
|20i3). We consider different approximations of Lagrangian pertur- 
bation theory: 1 st/Zeldovich, 2nd and curl-free 3rd order L PT We 
also use t he recent v e rsion of the spherical model (SC) (Bern ardeaul 
Il992h by iNevrinckl tOUh . Finally, we apply two improvements 
of the 2LPT model with *a(<?, z) = *2Lpt(<?, z). In the first 
one we include the ideal collapse limit in the 2LPT approxima- 
tion: ^b(9, z) = ^2LPTc(9, z). In the second one we addition- 
ally include the correction towards low densities using the spherical 
collapse model ^b(^, z) = ^sc(^, z). The relationship between 
the divergence of the displacement field and the linear overdensity 
field in various approximations is shown in Fig. [T] The linear be- 
haviour of the Zeldovich approximation is represented by the black 
line with slope of 1. The quadratic and the cubic behaviours of 
2nd and 3rd order LPT can be clearly seen in the scatter plots, left 
and middle panels, respectively. These were computed based on 
the 10th sample of our set of simulations. The SC model repre- 
sented by the green curve serves as a good proxy f or the mean rela - 
tion corresponding to the A^-body simulation (see lNevrinckll20 1 3h . 
The right panel shows the combined model 2LPT-SC which has 
been computed by taking the divergence of the resulting displace- 
ment field from Eq.|2] This relation has several characteristics. The 
mean fits well the SC model towards low overdensities. We note 
that the SC underestimates the overdensities at low values of the 
stretching parameter. Our model reproduces the results of the N - 
body simulation better in this regime (see Fig. 6 in lNevrinckll2013h . 
The stretching parameter truncates at stretching values of about 3-4 
simulating the collapse and suppressing shell crossing. The 2LPT- 
2LPTC model is similar to the 2LPT-SC, however, retaining the 
quadratic behaviour of the 2LPT model towards underdense re- 
gions. We move 256^ particles forward in time using the differ- 
ent displacement fields presented here. Fig. |2] shows the results of 
the 2LPT models, classic and improved ones, as compared to the 
A^-body simulation for the 10th sample. One can clearly see the 
strong shell crossing of 2LPT and how this is suppressed in knots 
and thick filaments in both improved models. A careful inspection 
shows also that the low density regions are better captured by the 
2LPT-SC model than by the 2LPT-2LPTC overdensity fields. It 
is remarkable how even small filaments which are present in the 
A/^-body simulation are also traced in the 2LPT-SC one. Another 
important characteristic of the improved 2LPT models is that they 
preserve power on large scales as opposed to the SC one (Fig. |3] 
shows the mean power- spectra of the 25 realisations for each ap- 
proximation). The improved 2LPT model has more power towards 
small scales than 2LPT showing that nonlinear structure formation 
is better captured. 

Finally, we assess the quality of our approach by computing 
the cross-power spectra between the various approximations and 
the A^-body solution. We find that our approach yields a higher cor- 
relation on all scales by choosing the appropriate threshold scale 
(see Fig. We perform a parameter study to determine the op- 
timal range of smoothing radii and find that the maximum cor- 
relation is achieved with rs = 4 — 5 Mpc (see right panel 
in Fig. |4]). This result is in agreement with previous works which 
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Figure 3. Mean power-spectra for different structure formation models us- 
ing the same set of initial conditions, black dashed: A^-body simulation; 
blue: LPT/Zeldovich; red: 2LPT; gray: 3LPT; green: SC; black: combined 
2LPT-SC approximation. 



also found th at 2LPT performs best down to scales of about 4 — 5 
M pc rsee lKitaura & Anguldl20ll : lKitaura et alJ2012l : lKitaural 
l2012h . When taking a smoothing scale of 6 /i ^ Mpc we start to get 
slightly worse results than 2LPT on very large scales k < OA 
Mpc. We note that the 2LPT-2LPTC model yields similar cross- 
power spectra to the 2LPT-SC model, indicating that shell crossing 
in high density regions reduces more dramatically the accuracy of 
2LPT. We have checked that the 2LPT approximation with a col- 
lapse truncation does not improve the results of the SC model. In 
fact, including the collapse limit significantly degrades the cross- 
correlation with respect to the classical 2LPT model. 

Furthermore, our results show that 3LPT does not improve the 
Zeldovich solution. We also find that the Zeldovich approximation 
performs better than 2LPT already at about k = 0.5 — 0.7/i Mpc~^ 
depending on the realisation. 



5 CONCLUSIONS 

We have presented a new approach to model structure formation 
augmenting Lagrangian perturbation theory (ALPT). Our method 
is based on splitting the displacement field into a long range and 
a short range component, given by 2LPT and the spherical col- 
lapse approximations, respectively. These are the advantages of our 
method: 

• It fits very well the relation between the divergence of the dis- 
placement field and the linear overdensity field. 

• It considerably enhances the cross-correlation on all scales be- 
tween the A^-body solution and previous approximations like the 
1st (Zeldovich), 2nd, 3rd order Lagrangian perturbation theory and 
the spherical model approximation. 

• It requires only one additional parameter which defines the 
transition from the long range to the short range regime. We find 
an optimal choice of this parameter to be in the range 4-5 Mpc 
for Gaussian smoothing. 

• Similar to 2LPT it preserves the power on large scales. 

• Since it is a one- step solver it is extremely fast and efficient as 
compared to A^-body simulations. 
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Figure 4. Normalised cross-power spectra XP(k) = {\Srec{k)6N\)ody{k)\)/{^yPrec{k)^/PNhody{k)) between the matter overdensity from an A^-body 
simulation and the corresponding simulation in Ist/Zeldovich (blue), 2nd (red), 3rd (gray) LPT, SC (green) and the combined 2LPT-SC (black) model. On 
the left: The mean cross-power spectra are given by the solid lines, the dashed line represents the 10th realisation of our set of simulations, the 1 sigma 
contours are given by the corresponding light shaded regions. The 2LPT-SC model was computed with a threshhold scale of rs =4 Mpc/h. On the right: 
same as left panel including the mean cross power- spectra for the following smoothing radii: rs = 3, 4, 5, 6 Mpc/h (dotted, solid, dashed-dotted and dashed, 
respectively). 



The range of applications is very wide. Due to its improved mod- 
eling of structure formation this approach could be used for the 
following purposes: 

• Set-up initial conditions for A/^-body sirn ulations (see 
ICrocce et alJl2006l : |jenkinsll201Ql : iKitaura et alJl2012h . 

• IVlake mock galaxy catalogues (see 2LPT based approaches 
IScoccimarro & Sheth 2002':' Manera et aLl20 12V which could be an 
improved tool for dark energy studies and cosmological parameter 
studies for present and future large-scale structure surveys. 

• Fast modeling of the Universe at high redshift to perform sta- 
tistical analysis (e.g. Lymana forest, 21 cm line) (s ee 21CMF AST 
which uses the Zeldovich approximation, Mesingeretal ]|201lh . 

• To analyse the cosmic structu re (cosmic web classification) 
and perform environmental studies dNuza et al.lboTi) . 

• For the inference of the primordial density fluctuations 
to perform baryon acoustic oscillation (BAO) reco nstructions 
( Eisenstein et al.l l2007h or constrained simulations ( Klypin et ajj 
l2003h . The method presented here can be directly implemented in 
KIGEN (lKitaural2012l : lKitaura et alj2012l) to re cover the initial c on- 
ditions to perform constrained simulations see iHeB etal.l(l2012h . 

We will address each of these issues in forthcoming publications. 
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